T=4;
ts=0.01;
N=T/ts;
%两个线性调频信号相加
n1=0:0.01:T/2-0.01;
x1=sin(2*pi*(0.1+10*n1).*n1);
x2=sin(2*pi*(0.1+5*n1).*n1);

my_result=mywvd(x1,x2,N/2);
mesh(my_result);

function myfft=mywvd(x,y,N)
x=hilbert(x);
dft_n=(N-1);
kn=zeros(dft_n,N);
for n=1:N
    if n<=N/2
        kmax=n-1;
    else
        kmax=N-n;
    end
    k_range=-kmax:kmax;
    indices=rem(N+k_range,N)+1;
    %第一种方法
%     for k = -kmax:kmax  
%         aa=N/2+k;
%         kn(aa,n)=x(n+k)*conj(y(n-k));
%     end
    %另一种方法:直接将序列生成标准形式，负的都位移到后面。
    kn(indices,n)=x(n+k_range).*conj(x(n-k_range));
end
myfft=fft(kn);
%用第一种方法需要加上下面循环语句
% for k = 1:dft_n
%     myfft(k,:)=myfft(k,:).*exp(-1i*pi*k);%移位回来。。。
% end
myfft=real(myfft);
%经过比较发现，第一种方法有一点不足，但是没有发现在哪里：负的值太多了
end
